#include "BWT.cu"
#include "cuda_header.h"
#include "mian_util.h"
#include <thrust/scan.h>
#include <thrust/reduce.h>
#include "thrust/sort.h"
#include "thrust/device_vector.h"
#include "thrust/host_vector.h"
#include <cuda.h>
//#include"OccValue_Lam.cu"

__device__ int BackwardSearch1Error3(uint3* result,uint4* range, unsigned int*state, int maxNumOutput,
		BWT * bwt, unsigned int* theQuery,
		unsigned int queryId,
		unsigned int lower, unsigned int upper, unsigned int mismatchCount, 
		unsigned int * sa_start, unsigned int * sa_end, 
		long* numAlign = NULL)
{
	unsigned int temp_start=(*sa_start),temp_end=*sa_end;
	int i;
	char c, misBase;
	int numResult = 0;

	for (i=upper-1; i >=(int)lower && *sa_start<=*sa_end; i--)
	{
		c = getCharFromQueries(theQuery, i);

		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{		
			if (misBase == c)
				continue;
			temp_start = *sa_start;
			temp_end = *sa_end;

			if(!bwtCalc(bwt,misBase,&temp_start,&temp_end))
				continue;

			//add pausedQuery to buffer	
			if(BackwardSearch(bwt, theQuery, 0,i, &temp_start, &temp_end))
			{
				result[numResult].x = queryId;
				result[numResult].y = temp_start;
				result[numResult].z = temp_end;
				numResult++;
				if(numResult + 1 > maxNumOutput) {
					(*state) |= ((i<<8)|misBase)&0xFFFF;
					return -1;
				}
			}
		}
		bwtCalc(bwt, c,sa_start,sa_end);
	}

	return numResult;
}



__device__ int ForwardSearch1Error3(uint3* result, uint4 * range, unsigned int* state, int maxNumOutput,
		BWT * revbwt, unsigned int * query, unsigned int queryId,
		unsigned int lower, unsigned int upper, unsigned int mismatchCount, 
		unsigned int * sa_rev_start, unsigned int * sa_rev_end, 
		unsigned int * sa_start, unsigned int * sa_end, 
		long* numAlign = NULL)
{
	unsigned int temp_start=*sa_start,temp_end=*sa_end, temp_rev_start=*sa_rev_start, temp_rev_end=*sa_rev_end;
	unsigned int occCount[ALPHABET_SIZE];
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	unsigned int i;
	char c, misBase;
	int numResult = 0;


	for (i=lower; i<upper && *sa_start<= *sa_end; i++) {

		bwtCalcAll(revbwt,sa_rev_start, sa_rev_end,occCount_start, occCount_end, occCount);
		c = getCharFromQueries(query, i);
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;

			temp_rev_start = *sa_rev_start;
			temp_rev_end = *sa_rev_end;

			bwtCalcWithOccValue(revbwt->countValue, misBase,occCount_start[misBase], occCount_end[misBase], 
					&temp_rev_start, &temp_rev_end);

			temp_end = *sa_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if (temp_start > temp_end)
				continue;


			if(ForwardSearch(revbwt, query, i+1,upper, &temp_rev_start, &temp_rev_end, &temp_start,&temp_end))
			{
				if(numResult + 1 <= maxNumOutput) {

					result[numResult].x = queryId;
					result[numResult].y = temp_start;
					result[numResult].z = temp_end;
					numResult++;
					if( numResult + 1 > maxNumOutput) {
						(*state) |= ((i<<8)|misBase)&0xFFFF;
						return -1;
					}
				}

				//ReportResult(mismatchCount+1, temp_start, temp_end, saInterval);
			}

		}

		bwtCalcWithOccValue(revbwt->countValue, c,occCount_start[c], occCount_end[c], sa_rev_start, sa_rev_end);

		*sa_end = *sa_end - occCount[c];
		*sa_start = *sa_end - (*sa_rev_end - *sa_rev_start);
	}

	return numResult;
}
/*
   __device__ int BWT1MismatchSearch_Backward(uint3* result, int maxNumOutput,
   unsigned int * query, unsigned int queryId,
   BWT *bwt, BWT *revbwt, long * numAlign = NULL)
   {
   unsigned int start= 1, end = bwt->textLength-1;
   char c;
   unsigned int i;
   int tempN = 0, numOutput=0;
   unsigned long long packedPattern = 0;
   int queryLen = QUERY_LENGTH;

   for (i = 0; i <bwt->lookupTable->tableSize ; i++) {
   packedPattern<<=2;
   c = getCharFromQueries(query, queryLen - TABLE_SIZE +i);
   packedPattern |= (c & 3);
   }
   start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
   end = bwt->lookupTable->table[packedPattern];

//      BackwardSearch(bwt, Queries, queryId, sizeX, QUERY_LENGTH-1, &start, &end);
if(start <= end && BackwardSearch(bwt, query, (QUERY_LENGTH / 2), queryLen-TABLE_SIZE, &start, &end))
{
tempN = BackwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
bwt, query, queryId,0,(QUERY_LENGTH / 2),0,
&start, &end, numAlign);

if( numAlign!=NULL &&  *numAlign == -1)
return 0;
if(tempN == -1) {
return -1;
} else {
numOutput += tempN;
}
}
return numOutput;
}


__device__ int BWT1MismatchSearch_Forward(uint3* result, int maxNumOutput,
unsigned int * query, unsigned int queryId,
BWT *bwt, BWT *revbwt, long * numAlign = NULL)
{

unsigned int start= 1, end = bwt->textLength-1, revStart, revEnd;
char c;
unsigned int i;

unsigned long long l_packedPattern = 0;
unsigned long long r_packedPattern = 0;
unsigned long long rev_packedPattern = 0;
int numOutput = 0;
for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
l_packedPattern<<=2;
c = getCharFromQueries(query, i);
l_packedPattern |= ( c & 3 );
}
r_packedPattern = l_packedPattern;
l_packedPattern <<= (TABLE_SIZE - revbwt->lookupTable->tableSize)*2;
r_packedPattern <<= (TABLE_SIZE -revbwt->lookupTable->tableSize)*2;
r_packedPattern |= (1<<(TABLE_SIZE - revbwt->lookupTable->tableSize)*2) - 1;

for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
rev_packedPattern<<=2;
c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
rev_packedPattern |= (c & 3);
}

start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
end = bwt->lookupTable->table[r_packedPattern];
revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
revEnd = revbwt->lookupTable->table[rev_packedPattern];

if (end - start < revEnd - revStart) {
	revStart =revEnd - (end - start );
} else if (end - start > revEnd - revStart) {
	start = end - (revEnd - revStart);
}

ForwardSearch(revbwt, query, revbwt->lookupTable->tableSize, (QUERY_LENGTH / 2), &revStart, &revEnd, &start, &end);

if(start <= end){
	int tempN = ForwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
			revbwt, query, queryId, (QUERY_LENGTH/2),QUERY_LENGTH,0,
			&revStart,&revEnd,&start, &end, numAlign);
	if(numAlign != NULL && (*numAlign) == -1)
		return 0;
	if(tempN == -1) {
		return -1;
	} else {
		numOutput += tempN;
	}
}
return numOutput;
}

__global__ void BWT1MismatchBackward_Kernel(uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * Queries, unsigned int queryNum, BWT *bwt, BWT *revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT * ptr_bwt = bwt;
	BWT * ptr_revbwt = revbwt;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		numOutputCurQuery = 0;

		numOutputCurCase = BWT1MismatchSearch_Backward(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				Queries + QUERY_SIZE_IN_WORD*idx, idx, ptr_bwt, ptr_revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;
	}
}


__global__ void BWT1MismatchForward_Kernel(uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * Queries, unsigned int queryNum, BWT *bwt, BWT *revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT * ptr_bwt = bwt;
	BWT * ptr_revbwt = revbwt;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		numOutputCurQuery = 0;

		numOutputCurCase = BWT1MismatchSearch_Forward(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				Queries + QUERY_SIZE_IN_WORD*idx, idx, ptr_bwt, ptr_revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;
	}
}
	*/
__device__ int BWT2MismatchSearch_CaseA3(uint3* result, uint4* range, unsigned int* state, int maxNumOutput,
		unsigned int * query, unsigned int queryId, 
		BWT *bwt, BWT *revbwt)
{
	unsigned int temp_start=0,temp_end=0;
	unsigned long long packedPattern = 0;

	unsigned int i;
	char c, misBase;
	unsigned int start, end;
	for (i = 0; i < TABLE_SIZE ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(query, QUERY_LENGTH - TABLE_SIZE + i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];


	int numOutput = 0;
	int tempN = 0;

	if(start <= end &&  BackwardSearch(bwt, query, sizeY, QUERY_LENGTH - TABLE_SIZE, &start, &end))
	{
		for (i=sizeY-1; i>=1 && start<=end; i--)
		{
			c = getCharFromQueries(query, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
			{
				if( misBase == c)
					continue;
				temp_start = start;
				temp_end = end;
				if(bwtCalc(bwt,misBase,&temp_start,&temp_end)){
					tempN = BackwardSearch1Error3(result + numOutput, range, state, maxNumOutput - numOutput,
							bwt, query, queryId,0,i,1, &temp_start, &temp_end);
					if(tempN == -1) {
						(*state) <<= 16;
						(*state) |= (((i<<8)|misBase)&0xFFFF);
						range->x = temp_start;
						range->y = temp_end;
						range->z = 0;
						range->w = 0;
						return -1;
					} else {
						numOutput += tempN;
					}
				}
			}
			bwtCalc(bwt, c, &start, &end);
		}

	}

	return numOutput;
}

__device__ int BWT2MismatchSearch_CaseA_Continue(uint3* result, uint4* range, unsigned int* state, int maxNumOutput, 
		unsigned int * query, unsigned int queryId,  
		BWT *bwt, BWT *revbwt) 
{ 
	unsigned int temp_start=0,temp_end=0; 
	unsigned long long packedPattern = 0; 

	unsigned int i; 
	char c, misBase, c0; 
	unsigned int start, end; 

	start = range->x;
	end = range->y; 

	c0 = (*state) & 0x3;
	i = ((*state) >> 8)0xFF; 
	int numOutput = 0; 
	int tempN = 0; 

	c = getCharFromQueries(query, i); 
	for(misBase = 0; misBase < ALPHABET_SIZE; misBase++) {
		if(misBase <= c0 || misBase == c)
			continue;
		temp_start = start; 
		temp_end = end; 
		if(bwtCalc(bwt,misBase,&temp_start,&temp_end)){ 
			tempN = BackwardSearch1Error3(result + numOutput, range, state, maxNumOutput - numOutput, 
					bwt, query, queryId,0,i,1, &temp_start, &temp_end); 
		}
	}	 
	for (i=i-1; i>=1 && start<=end; i--) 
	{ 
		c = getCharFromQueries(query, i); 
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++) 
		{ 
			if( misBase == c) 
				continue; 
			temp_start = start; 
			temp_end = end; 
			if(bwtCalc(bwt,misBase,&temp_start,&temp_end)){ 
				tempN = BackwardSearch1Error3(result + numOutput, range, state, maxNumOutput - numOutput, 
						bwt, query, queryId,0,i,1, &temp_start, &temp_end); 
				if(tempN == -1) { 
					(*state) <<= 16; 
					(*state) |= (((i<<8)|misBase)&0xFFFF); 
					range->x = temp_start; 
					range->y = temp_end; 
					range->z = 0; 
					range->w = 0; 
					return -1; 
				} else { 
					numOutput += tempN; 
				} 
			} 
		} 
		bwtCalc(bwt, c, &start, &end); 
	} 

	return numOutput; 
} 

__device__ int BWT2MismatchSearch_CaseB3(uint3* result, uint4*range, unsigned int *state,int maxNumOutput,
		unsigned int * query, unsigned int queryId, 
		BWT *bwt, BWT *revbwt)
{
	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;

	unsigned int occCount[ALPHABET_SIZE];
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	unsigned long long l_packedPattern = 0;
	unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;
	int numOutput = 0;

	l_packedPattern = 0;
	r_packedPattern = 0;
	rev_packedPattern = 0;
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		l_packedPattern<<=2;
		c = getCharFromQueries(query, i);
		l_packedPattern |= (c & 3);
	}
	r_packedPattern = l_packedPattern;
	l_packedPattern <<= (bwt->lookupTable->tableSize -  revbwt->lookupTable->tableSize)*2;
	r_packedPattern <<= (bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2;
	r_packedPattern |= (1<<(bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2) - 1;

	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
	end  = bwt->lookupTable->table[r_packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	if (end - start<revEnd-revStart)
	{
		revStart = revEnd - (end - start);
	} else if (end - start > revEnd - revStart)
	{
		start = end - (revEnd - revStart);
	}

	if (start <= end && ForwardSearch(revbwt, query, revbwt->lookupTable->tableSize,sizeY,
				&revStart,&revEnd,&start,&end))
	{
		for (i=sizeY; i<sizeZ && start<=end; i++)
		{
			bwtCalcAll(revbwt, &revStart, &revEnd, occCount_start, occCount_end, occCount);
			c = getCharFromQueries(query, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
			{
				if(misBase == c)
					continue;
				temp_rev_start = revStart;
				temp_rev_end = revEnd;

				bwtCalcWithOccValue(revbwt->countValue,misBase,
						occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

				temp_end = end - occCount[misBase];
				temp_start = temp_end - (temp_rev_end - temp_rev_start);

				if(temp_start <= temp_end){

					int tempN = ForwardSearch1Error3(result + numOutput, range, state, maxNumOutput - numOutput,
							revbwt, query, queryId, i+1,sizeZ, 1, 
							&temp_rev_start, &temp_rev_end,&temp_start, &temp_end);
					if(tempN == -1) {
						range.x = temp_start;
						range.y = temp_end;
						range.z = temp_rev_start;
						range.w = temp_rev_end;
						(*state) <<= 16;
                                                (*state) |= (((i<<8)|misBase)&0xFFFF);
						return -1;
					}
					else {
						numOutput += tempN;
					}

				}
			}

			bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

			end = end - occCount[c];
			start = end - (revEnd - revStart);
		}
	}

	return numOutput;
}

__device__ int BWT2MismatchSearch_CaseB_Continue(uint3* result, uint4*range, unsigned int *state,int maxNumOutput, 
                unsigned int * query, unsigned int queryId,  
                BWT *bwt, BWT *revbwt) 
{ 
        unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0; 
 
        unsigned int occCount[ALPHABET_SIZE]; 
        unsigned int occCount_start[ALPHABET_SIZE]; 
        unsigned int occCount_end[ALPHABET_SIZE]; 
 
 
        unsigned int i,pos0; 
        char c, misBase, c0; 
        unsigned int start, end, revStart, revEnd; 
        int numOutput = 0; 
 
	start = range->x;
	end = range->y;
	revStart = range->z;
	revEnd = range->w;

	c0 = (*state) & 0x3;
	pos0 = ((*state) >> 8) & 0xFF; 
	
        c = getCharFromQueries(query, i); 
             for (i=pos0; i<sizeZ && start<=end; i++) 
                { 
                        bwtCalcAll(revbwt, &revStart, &revEnd, occCount_start, occCount_end, occCount); 
                        c = getCharFromQueries(query, i); 
                        for(misBase=0; misBase<ALPHABET_SIZE; misBase++) 
                        { 
				if(i==pos0 && misBase <= c0)
					continue;
                                if(misBase == c) 
                                        continue; 
                                temp_rev_start = revStart; 
                                temp_rev_end = revEnd; 
 
                                bwtCalcWithOccValue(revbwt->countValue,misBase, 
                                                occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end); 
 
                                temp_end = end - occCount[misBase]; 
                                temp_start = temp_end - (temp_rev_end - temp_rev_start); 
 
                                if(temp_start <= temp_end){ 
 
                                        int tempN = ForwardSearch1Error3(result + numOutput, range, state, maxNumOutput - numOutput, 
                                                        revbwt, query, queryId, i+1,sizeZ, 1,  
                                                        &temp_rev_start, &temp_rev_end,&temp_start, &temp_end); 
                                        if(tempN == -1) { 
                                                range.x = temp_start; 
                                                range.y = temp_end; 
                                                range.z = temp_rev_start; 
                                                range.w = temp_rev_end; 
						(*state) <<= 16;
                                                (*state) |= (((i<<8)|misBase)&0xFFFF);
                                                return -1; 
                                        } 
                                        else { 
                                                numOutput += tempN; 
                                        } 
 
                                } 
                        } 
 
                        bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd); 
 
                        end = end - occCount[c]; 
                        start = end - (revEnd - revStart); 
                } 
 
        return numOutput; 
} 


__device__ int BWT2MismatchSearch_CaseC3(uint3* result,uint4* range, unsigned int* state, int maxNumOutput,
		unsigned int * query, unsigned int queryId, 
		BWT *bwt, BWT *revbwt)
{

	int numOutput = 0;

	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;

	unsigned int occCount[ALPHABET_SIZE];
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	//	unsigned int sizeX = floor((double)QUERY_LENGTH * 0.3f);
	//	unsigned int sizeY = sizeX * 2;
	//	unsigned int sizeZ = QUERY_LENGTH;

	unsigned long long packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	unsigned long long ALLONE = (1<<(bwt->lookupTable->tableSize*2))-1;
	unsigned long long mask;

	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;

	packedPattern = 0;
	rev_packedPattern = 0;

	for (i = 0; i <bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(query, i);
		packedPattern |= (c & 3);
	}
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3 );
	}

	//for error happens in lookup range
	for(i=sizeX; i<bwt->lookupTable->tableSize; i++)
	{
		c = getCharFromQueries(query, i);
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{
			if(misBase == c)
				continue;
			unsigned long long err_packedPattern = packedPattern;
			unsigned long long err_rev_packedPattern = rev_packedPattern;

			unsigned int bitPos = (bwt->lookupTable->tableSize-i-1)*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_packedPattern |= (3 << bitPos);
			err_packedPattern &= mask;

			bitPos = i*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_rev_packedPattern |= (3 << bitPos);
			err_rev_packedPattern &= mask;

			start = err_packedPattern ? bwt->lookupTable->table[err_packedPattern-1]+1 : 1;
			end = bwt->lookupTable->table[err_packedPattern];
			revStart = err_rev_packedPattern ? revbwt->lookupTable->table[err_rev_packedPattern-1]+1 : 1;
			revEnd = revbwt->lookupTable->table[err_rev_packedPattern];

			//Lookup SA correction...
			if (end - start < revEnd - revStart) {
				revStart = revEnd - ( end - start);
			}else if ( end - start > revEnd - revStart) {
				revStart = revEnd -( revEnd - revStart);
			}

			if(start > end)
				continue;
			unsigned int pos = revbwt->lookupTable->tableSize;
			ForwardSearch(revbwt, query, pos, sizeY, &revStart, &revEnd, &start, &end);

			if( start <= end ){

				int tempN = ForwardSearch1Error3(result + numOutput, range, state, maxNumOutput - numOutput,
						revbwt,  query, queryId,sizeY, sizeZ, 1, 
						&revStart,&revEnd,&start, &end);

				if(tempN == -1) {
					range.x = start;
					range.y = end;
					range.z = revStart;
					range.w = revEnd;
						(*state) <<= 16;
                                                (*state) |= (((i<<8)|misBase)&0xFFFF);
					return -1;
				} else {
					numOutput += tempN;
				}
			}
		}
	}




	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	revStart = packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];





	//Lookup SA correction...
	if (end - start < revEnd - revStart) {
		revStart = revEnd - ( end - start);
	} else if (end - start > revEnd - revStart) {
		start = end - ( revEnd - revStart);
	}

	if(sizeX > revbwt->lookupTable->tableSize){
		ForwardSearch(revbwt, query, TABLE_SIZE ,sizeX,
				&revStart,&revEnd,&start,&end);
		i = sizeX;
	}
	else 
		i = TABLE_SIZE;
	for (; i<sizeY && start<=end ; i++)
	{
		bwtCalcAll(revbwt, &revStart, &revEnd,occCount_start, occCount_end, occCount);
		c = getCharFromQueries(query, i);
		for(misBase=0; misBase < ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = start;
			temp_end = end;
			temp_rev_start = revStart;
			temp_rev_end = revEnd;

			bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

			temp_end = temp_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if(temp_start > temp_end)
				continue;
			if(ForwardSearch(revbwt, query, i+1,sizeY,&temp_rev_start,&temp_rev_end,&temp_start,&temp_end) 
					&& temp_start <= temp_end) {

				int tempN = ForwardSearch1Error3(result + numOutput, range, state, maxNumOutput - numOutput,
						revbwt, query, queryId,sizeY,sizeZ,1,&temp_rev_start,
						&temp_rev_end,&temp_start, &temp_end);

				if(tempN == -1) {
					range.x = temp_start;
					range.y = temp_end;
					range.z = temp_rev_start;
					range.w = temp_rev_end;
						(*state) <<= 16;
                                                (*state) |= (((i<<8)|misBase)&0xFFFF);
					return -1;
				} else {
					numOutput += tempN;
				}
			}
		}

		bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

		end = end - occCount[c];
		start = end - (revEnd - revStart);
	}

	return numOutput;
}




__device__ int BWT2MismatchSearch_CaseD3(uint3* result, int maxNumOutput,uint4* range, unsigned int *state,
		unsigned int * Queries, unsigned int queryId, BWT *bwt, BWT *revbwt)
{

	int numOutput = 0;

	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
	unsigned int occCount[ALPHABET_SIZE];
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	//	unsigned int sizeX = floor((double)QUERY_LENGTH * 0.3f);
	//	unsigned int sizeY = sizeX * 2;
	//	unsigned int sizeZ = QUERY_LENGTH;

	unsigned long long packedPattern = 0;
	//unsigned long long l_packedPattern = 0;
	//unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;
	//unsigned long long rev_l_packedPattern = 0;
	//unsigned long long rev_r_packedPattern = 0;

	unsigned long long ALLONE = (1<<(bwt->lookupTable->tableSize*2))-1;
	unsigned long long mask;

	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;

	packedPattern = 0;
	//l_packedPattern = 0;
	//r_packedPattern = 0;
	rev_packedPattern = 0;
	//rev_l_packedPattern = 0;
	//rev_r_packedPattern = 0;

	for (i = 0; i < bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(Queries, sizeX+i);
		packedPattern |= (c & 3);
	}
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(Queries, sizeX+revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	//For error happen in lookup range....
	for (i = sizeY - sizeX; i< bwt->lookupTable->tableSize ;i++) {
		c = getCharFromQueries(Queries, sizeX + i);
		for (misBase=0;misBase<4;misBase++) {
			if (misBase == c) continue;
			unsigned long long err_packedPattern = packedPattern;
			unsigned long long err_rev_packedPattern = rev_packedPattern;
			unsigned int bitPos = (bwt->lookupTable->tableSize-i-1)*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_packedPattern |= (3 << bitPos);
			err_packedPattern &= mask;

			bitPos = i*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_rev_packedPattern |= (3 << bitPos);
			err_rev_packedPattern &= mask;

			start = err_packedPattern ? bwt->lookupTable->table[err_packedPattern-1]+1 : 1;
			end = bwt->lookupTable->table[err_packedPattern];
			revStart = err_rev_packedPattern ? revbwt->lookupTable->table[err_rev_packedPattern-1]+1 : 1;
			revEnd = revbwt->lookupTable->table[err_rev_packedPattern];

			//Lookup SA correction...
			if (end - start < revEnd - revStart) {
				revStart = revEnd - ( end - start );
			} else if ( end - start >revEnd - revStart) {
				start = end - (revEnd - revStart);
			}
			if( start <= end && ForwardSearch(revbwt, Queries, sizeX+revbwt->lookupTable->tableSize, QUERY_LENGTH, 
						&revStart, &revEnd, &start, &end )) {


				int tempN = BackwardSearch1Error3(result + numOutput,range, state, maxNumOutput - numOutput,
						bwt, Queries, queryId, 0, sizeX, 1, &start,&end);
				if(tempN == -1) {
					range.x = start;
					range.y = end;
					range.z = revStart;
					range.w = revEnd;
					(*state) <<= 16;
					(*state) |= (((i<<8)|misBase)&0xFFFF);
					return -1;
				} else {
					numOutput += tempN;
				}

				//BackwardSearch1Error(bwt, Queries, queryId, 0, sizeX, 1, &start,&end,saInterval);
			}
		}
	}

	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end-start<revEnd-revStart) {
		revStart=revEnd-(end-start);
	} else if (end-start>revEnd-revStart) {
		start=end-(revEnd-revStart);
	}

	if(sizeY - sizeX > TABLE_SIZE) {
		ForwardSearch(revbwt, Queries, sizeX + TABLE_SIZE ,sizeY,
				&revStart,&revEnd,&start,&end);
		i = sizeY;

	}
	else
		i = sizeX + TABLE_SIZE;
	for (; i<sizeZ && start<=end; i++)
	{
		c = getCharFromQueries(Queries, i);
		bwtCalcAll(revbwt, &revStart, &revEnd,occCount_start, occCount_end, occCount);
		for (misBase = 0; misBase<ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = start;
			temp_end = end;
			temp_rev_start = revStart;
			temp_rev_end = revEnd;

			bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

			temp_end = temp_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if (temp_start <= temp_end)
				ForwardSearch(revbwt, Queries, i+1,sizeZ, &temp_rev_start, &temp_rev_end, &temp_start, &temp_end);
			if(temp_start <= temp_end){

				int tempN = BackwardSearch1Error3(result + numOutput,range, state, maxNumOutput - numOutput,
						bwt, Queries, queryId,0,sizeX,1,&temp_start, &temp_end);
				if(tempN == -1) {
					range.x = temp_start;
					range.y = temp_end;
					range.z = temp_rev_start;
					range.w = temp_rev_end;
					(*state) <<= 16;
					(*state) |= (((i<<8)|misBase)&0xFFFF);		
					return -1;
				} else {
					numOutput += tempN;
				}

				//BackwardSearch1Error(bwt, Queries, queryId,0,sizeX,1,&temp_start, &temp_end, saInterval);
			}

		}
		bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

		end = end - occCount[c];
		start = end - (revEnd - revStart);
	}


	return numOutput;
}





__global__ void BWT2MismatchSearchCaseA_continue_kernel(	
		uint3* d_result,uint4 *d_range, unsigned int *d_state, int* d_numOutput, 
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryA[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/

		if(d_numOutput[idx] > 0)
			continue;
		numOutputCurQuery = 0;	
		numOutputCurCase = BWT2MismatchSearch_CaseA3(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				d_range + idx, d_state + idx,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], 
				&bwt, &revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;		

	}
}



__global__ void BWT2MismatchSearchCaseB_continue_kernel(	
		uint3* d_result, uint4 *d_range, unsigned int* d_state,
		int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	//extern __shared__ unsigned int s_queryB[];
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryB[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/
		if(d_numOutput[idx] > 0)
			continue;

		numOutputCurQuery = 0;	
		numOutputCurCase = BWT2MismatchSearch_CaseB3(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				d_range + idx, d_state + idx, 
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;		
	}
}



__global__ void BWT2MismatchSearchCaseC_continue_kernel(	
		uint3* d_result, uint4 * d_range, unsigned int* d_state,
		int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	//extern __shared__ unsigned int s_queryC[];
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryC[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/
		if(d_numOutput[idx] > 0)
			continue;

		numOutputCurQuery = 0;

		numOutputCurCase = BWT2MismatchSearch_CaseC3(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				d_range + idx, d_state + idx,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;		
	}
}




__global__ void BWT2MismatchSearchCaseD_continue_kernel(	
		uint3* d_result, uint4 * d_range, unsigned int* d_state,
		int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	//extern __shared__ unsigned int s_queryD[];
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryD[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/
		if(d_numOutput[idx] > 0)
			continue;

		numOutputCurQuery = 0;
		numOutputCurCase = BWT2MismatchSearch_CaseD3(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				d_range + idx, d_state + idx,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;		
	}
}

__global__ void BWT2MismatchSearch_4Case_kernel(uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	//extern __shared__ unsigned int s_queryD[];
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryD[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/

		numOutputCurQuery = 0;
		numOutputCurCase = 0;
		numOutputCurCase = BWT2MismatchSearch_CaseA2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt);
		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}

		numOutputCurCase = BWT2MismatchSearch_CaseB2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt);
		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}

		numOutputCurCase = BWT2MismatchSearch_CaseC2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt);
		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}

		numOutputCurCase = BWT2MismatchSearch_CaseD2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;
	}
}

